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Abstract In this contribution we review recent efforts 
on investigations of the effect of (apparent) boundary 
slip by utilizing lattice Boltzmann simulations. We demon- 
strate the applicability of the method to treat funda- 
mental questions in microfluidics by investigating fluid 
flow in hydrophobic and rough microchannels as well as 
over surfaces covered by nano- or microscale gas bub- 
bles. 

Keywords Apparent and intrinsic slip • rough and 
hydrophobic surfaces • lattice Boltzmann simulations 
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1 Introduction 

During the last few decades the miniaturization of tech- 
nical devices down to submicrometric sizes has made 
considerable progress. In particular, so-called micro- 
electro-mechianical systems (MEMS) became available 
for chemical, biological and technical applications lead- 
ing to the rise of "microfluidics" about 20 years ago [1]. 
A wide variety of microfluidic systems including gas 
chromatography systems, electrophoretic separation sys- 
tems, micromixers, DNA amplifiers, and chemical reac- 

Jens Harting 

Department of Applied Physics, TU Eindhoven 

Den Dolech 2, 5600MB Eindhoven, The Netherlands 

Institute for Computational Physics, University of Stuttgart, 

Pfaffenwaldring 27, 70569 Stuttgart, Germany 

E-mail: j.harting@tuc.nl 

Christian Kunert 

Institute for Computational Physics, University of Stuttgart, 
Pfaffenwaldring 27, 70569 Stuttgart, Germany 
E-mail: kuni@icp . uni-stuttgart .de 

Jari Hyvaluoma 

Department of Physics, University of Jyvaskyla, FI-40014 

Jyvaskyla, Finland 

E-mail: jari.hyvaluoma@jyu.fi 



tors were developed. Next to those "practical applica- 
tions", microfluidics was used to answer fundamental 
questions in physics including the behavior of single 
molecules or particles in fluid flow or the validity of 
the no-slip boundary condition [1,2]. The latter is the 
focus of the current review and is investigated in detail 
by mesoscopic computer simulations. 

Reynolds numbers in microfluidic systems are usu- 
ally small, i.e., usually below 0.1. In addition, due to the 
small scales of the channels, the surface to volume ratio 
is high causing surface effects like wettability or surface 
charges to be more important than in macroscopic sys- 
tems. Also, the mean free path of a fluid molecule might 
be of the same order as the characteristic length scale 
of the system. For gas flows, this effect can be char- 
acterized by the so-called Knudsen number [3]. While 
the Knudsen number provides a good estimate for when 
to expect rarefaction effects in gas flows, for liquids one 
would naively assume that its velocity close to a surface 
always corresponds to the actual velocity of the surface 
itself. This assumption is called the no-slip boundary 
condition and can be counted as one of the generally ac- 
cepted fundamental concepts of fluid mechanics. How- 
ever, this concept was not always well accepted. Some 
centuries ago there were long debates about the veloc- 
ity of a Newtonian liquid close to a surface and the ac- 
ceptance of the no-slip boundary condition was mostly 
due to the fact that no experimental violations could 
be found, i.e., a so-called boundary slip could not be 
detected. 

In recent years, it became possible to perform very 
well controlled experiments that have shown a violation 
of the no-slip boundary condition in sub-micron sized 
geometries. Since then, mostly experimental [2, 4, 5, 5— 
10], but also theoretical works [11,12], as well as com- 
puter simulations [13-16, 31] have been performed to 
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improve our understanding of boundary slip. The topic 
is of fundamental interest because it has practical conse- 
quences in the physical and engineering sciences as well 
as for medical and industrial applications. Interestingly, 
also for gas flows, often a slip length much larger than 
expected from classical theory can be observed. Exten- 
sive reviews of the slip phenomenon have recently been 
published by Lauga et al. [2], Neto et al. [18], as well as 
Bocquet and Barrat [19]. 

The reason for such findings is that the behavior of 
a fluid close to a solid interface is very complex and 
involves the interplay of many physical and chemical 
properties. These include the wettability of the solid, 
the shear rate or flow velocity, the bulk pressure, the 
surface charge, the surface roughness, as well as im- 
purities and dissolved gas. Since all those quantities 
have to be determined very precisely, it is not sur- 
prising that our understanding of the phenomenon is 
still very unsatisfactory. Due to the large number of 
different parameters, a significant dispersion of the re- 
sults can be observed for almost similar systems [2, 18]. 
For example, observed slip lengths vary between a few 
nanometres [20] and micrometers [5] and while some 
authors find a dependence of the slip on the flow veloc- 
ity [4,7,21], others do not [5,6]. 

A boundary slip is typically quantified by the so- 
called slip length (3 - a concept that was already pro- 
posed by Navier in 1823. He introduced a boundary 
condition where the fluid velocity at a surface is propor- 
tional to the shear rate at the surface [22] (at x = x ), 
i.e., 



In other words, the slip length {3 can be defined as the 
distance from the surface where the relative flow veloc- 
ity vanishes. Assuming a typical Poiscuillc setup con- 
sisting of a pressure driven flow of an incompressible 
liquid between two infinite planes, the velocity in flow 
direction (y z ) at position x between the planes is given 
by 

v z (x) = ±^[d 2 -x 2 -2dP], (2) 

where 2d is the distance between the planes, and [i the 
dynamic viscosity. dP/dz is the pressure gradient. In 
contrast to a no-slip formulation, the last term in Eq.|2] 
linearly depends on the slip length (3. 

Most recent computer simulations apply molecular 
dynamics and report increasing slip with decreasing liq- 
uid density [23,24] or liquid-solid interactions [15,25], 
while slip decreases with increasing pressure [14] . These 
simulations are usually limited to a few tens of thou- 
sands of particles, length scales of a few nanometres 



and time scales of nanoseconds. Also, shear rates are 
usually some orders of magnitude higher than in any 
experiment [2]. Due to the small accessible time and 
length scales of molecular dynamics simulations, meso- 
scopic simulation methods such as the lattice Boltz- 
mann method arc well applicable for the simulation of 
microfluidic experiments. 

The experimental investigation of apparent slip can 
be based on different setups: either a fluid is pumped 
through a microchanncl and the measured mass flow 
rate at the end of the channel is compared to the theo- 
retical value with no slip boundary conditions. From the 
deviation of the two values, the magnitude of slip can 
be computed [26] . Another possibility is to measure the 
slip length directly using optical methods like particle 
image velocimetry (PIV). Very popular is the modifica- 
tion of an atomic force microscope (AFM) by adding a 
silicon sphere to the tip of the cantilever. While moving 
the sphere towards the boundary, the required force is 
measured. It is possible to measure the amount of slip 
at the wall by comparing the force needed to move the 
sphere with its theoretical value [10,27]. 

During the last few years, the substantial scientific 
research invested in the slip phenomenon has lead to a 
more clear picture which can be summarized as follows: 
one can argue that many surprising results published 
were only due to artefacts or misinterpretation of ex- 
periments. In general, there seems to be an agreement 
within the community that slip lengths larger than a 
few nanometers can usually be referred to as "apparent 
slip" and are often caused by experimental artefacts. 
Small slip lengths are experimentally even harder to 
determine and require sophisticated setups such as the 
modified atomic force microscopes as described above. 
Here, small variations of the apparatus such as choos- 
ing a different shape of the cantilever or modifying the 
control circuit of the sample holder can lead to substan- 
tial variation of the measurements. Also, the theoreti- 
cal equations correlating the measured force to the slip 
length are only valid for perfect surfaces and infinitely 
slow oscillations of the sphere. Therefore, it is of impor- 
tance to perform computer simulations which have the 
advantage that most parameters can be changed inde- 
pendently without modifying anything else. Thus, the 
influence of every single modification can be studied in 
order to present estimates of expected slip lengths. 

2 Apparent slip in hydrophobic microchannels 

The simulation method used to study microfluidic de- 
vices has to be chosen carefully. While Navier-Stokes 
solvers are able to cover most problems in fluid dynam- 
ics, they lack the possibility to include the influence 
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of molecular interactions as needed to model bound- 
ary slip. Molecular dynamics simulations (MD) are the 
best choice to simulate the fluid-wall interaction, but 
the computer power today is not sufficient to simu- 
late length and time scales necessary to achieve or- 
ders of magnitude which are relevant for experiments. 
However, boundary slip with a slip length (3 of the or- 
der of many molecular diameters a has been studied 
with molecular dynamics simulations by various au- 
thors [8,15,16,28,29]. 

The current contribution focuses on numerical in- 
vestigations of the slip phenomenon by means of lat- 
tice Boltzmann simulations. While an emphasis is put 
on reviewing our own contributions to the field, the 
achievements of other groups are commonly referred to. 
However, it should be noticed that while a large num- 
ber of groups utilizes the lattice Boltzmann technique 
to investigate microfluidic problems, only a very small 
number of researchers is actually applying the method 
to studying slippage. Even though interactions have to 
be described on a mesoscopic scale, this is surprising 
since mesoscopic simulation methods offer a closer re- 
lation to experimentally relevant time and length scales 
than microscopic techniques such as molecular dynam- 
ics. 

In the lattice Boltzmann method, one discretizes the 
Boltzmann kinetic equation 
d 

— + v\7 x rj(x,v,t) =n (3) 

on a lattice. The Boltzmann kinetic equation describes 
the evolution of the single particle probability density 
7/(x, v, t), where x is the position, v the velocity, and t 
the time. The derivatives represent simple propagation 
of a single particle in real and velocity space whereas 
the collision operator Q. takes into account molecular 
collisions in which a particle changes its momentum due 
to a collision with another particle. To represent the 
correct physics, the collision operator should conserve 
mass and momentum, and should be Galilei invariant. 
By performing a Chapman Enskog procedure, it can be 
shown that such a collision operator i~2 reproduces the 
Navier-Stokes equation [30]. In the lattice Boltzmann 
method the time t, the position x, and the velocity v 
are discretized. 

A few groups have applied the lattice Boltzmann 
method for the simulation of microflows and to study 
boundary slip. A popular approach is to introduce slip 
by generalizing the no-slip bounce back boundary con- 
ditions in order to allow specular reflections with a 
given probability [13,31-33], or to apply diffuse scat- 
tering [34-36]. It has been shown by Guo et al. that 
these approaches are virtually equivalent [37]. Another 
possibility is to modify the fluid's viscosity, i.e., the fluid 



viscosity is modified due to local density variations in 
order to model slip [38]. In both cases, the parame- 
ters determining the properties at the boundaries are 
"artificial" parameters and they do not have any ob- 
vious physical meaning. Therefore, they are not easily 
mappable to experimentally available values. We model 
the interaction between hydrophobic channel walls and 
the fluid by means of a multi-phase lattice Boltzmann 
model. Our approach overcomes this problem by apply- 
ing a mesoscopic force between the walls and the fluid. 
A similar approach is used by Zhu et al. [39], Benzi et 
al. [40], and Zhang et al. [41]. This force applied at the 
boundary can be linked to the contact angle which is 
commonly used by experimentalists to quantitatively 
describe the wettability of a material [42,43]. 

The simulation method and our implementation of 
boundary conditions are described as follows. A multi- 
phase lattice Boltzmann system can be represented by 
a set of equations 

r,f (x + c u t + 1) - <(x,i) = Qf , i = 0, 1, . . . ,b , (4) 

where ry"(x,t) is the single-particle distribution func- 
tion, indicating the amount of species a with velocity 
Ci, at site x on a D-dimensional lattice of coordination 
number b (D3Q19 in our implementation), at time-step 
t. This is a discretized version of equation ^ without 
external forces F for a number of species a. For the 
collision operator flf we choose the Bhatnagar-Gross- 
Krook (BGK) form [44] 

1 
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(5) 



where r Q is the mean collision time for component a 
and determines the kinematic viscosity 
2r Q — 1 

,° = ^V- (•) 

of the fluid. The relaxation time t q is kept constant at 
1.0 in this study. The system relaxes to an equilibrium 
distribution 77" e<3 ' which can be derived imposing re- 
strictions on the microscopic processes, such as explicit 
mass and momentum conservation for each species. In 
our implementation we choose for the equilibrium dis- 
tribution function 
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(7) 



which is a polynomial expansion of the Maxwell dis- 
tribution. Ci are the velocity vectors pointing to neigh- 
bouring lattice sites and Q are the lattice weights result- 
ing from the velocity space discretization. c s = 1/V3 is 
the speed of sound for the D3Q19 lattice. The macro- 
scopic values can be derived from the single-particle 
distribution function ?7"(x, t), i.e., the density 77" (x, t) 
of the species a at lattice site x is the sum over the 
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distribution functions r)"(x,t) for all lattice velocities 

Ci, 

rj a ( X ,t) = J2v?(^t). (8) 

i 

u Q (x, t) is the macroscopic velocity of the fluid, defined 
as 

r, Q (x,t)u Q (x,t)^^r,r(x,t)c l . (9) 

i 

Interactions between different fluid species are intro- 
duced following Shan and Chen as a mean field body 
force between nearest neighbors [45,46], 

F Q (x, t) ee -#*(x, t) £ g a& J2 ^( x '' *)( x ' - x ) . ( 10 ) 

a x' 

where ijj a (x,t) = (l_ e -'7 0, (*»*)/'?o) is the so-called effec- 
tive mass with r/o being a reference density that is set 
to 1 in our case [45]. g a& is a force coupling constant, 
whose magnitude controls the strength of the interac- 
tion between component a and a. The dynamical effect 
of the force is realized in the BGK collision operator ([5| 
by adding an increment <5u Q = r Q F Q /rj a to the velocity 
u in the equilibrium distribution function ^ . A repul- 
sive potential between surface and fluid can be used 
to model hydrophobic fluid-surface interactions. Such 
a potential is realized by attaching the imaginary fluid 
"density" jy wa11 to the first lattice site inside the wall. 
Only the distribution corresponding to the rest velocity 
is filled, while the remaining ones are kept at 0. As a 
result the only difference between r] wan and any other 
fluid packages on the lattice rj a is that the fluid corre- 
sponding to J7 wa11 is only taken into account for in the 
collision step and for the calculation of Eq. [10] but not 
in the propagation step. Therefore, we can adopt 7y wa11 
and the coupling constant <7a iWa ii in order to tune the 
fluid-wall interaction, g^waii is kept at 0.08 throughout 
this paper if not mentioned otherwise and all values are 
reported in lattice units. These parameters allow to sim- 
ulate a wide range of effective interactions without com- 
promising on numerical stability. Additionally, we apply 
second order correct mid-grid bounce back boundary 
conditions between the fluid and the surface which as- 
sures vanishing velocities at solid surfaces. Here, a dis- 
tribution function that would be advected into a solid 
node is simply reversed and advected into the opposite 
direction [30]. 

From molecular dynamics simulations it is known 
that the fluid-wall interactions causing a slip phenomenon 
usually take place within a few molecular layers of the 
liquid along the boundary surface [8,15,16,28]. Our 
coarse-grained fluid wall interaction acts on the length 
scale of one lattice constant and does not take the molec- 
ular details into account. Therefore, coarse-grained im- 
plementations based on the lattice Boltzmann method 



are only able to reproduce an averaged effect of the in- 
teraction and cannot fully resolve the correct flow pro- 
file very close to the wall and below the resolution of a 
single lattice spacing. However, in order to understand 
the influence of the hydrophobicity on experimentally 
observed apparent slip, it is fully sufficient to inves- 
tigate the flow behavior on more macroscopic scales 
as they are accessible for experimental investigation. 
Coarse-grained interaction models could be improved 
by a direct mapping of data obtained from MD simu- 
lations to the coupling constant g Q . W aii allowing a di- 
rect comparison of the influence of liquid-wall interac- 
tions on the detected slip [47]. Similar approaches are 
known from quantitative comparisons of lattice Boltz- 
mann and molecular dynamics simulations in the liter- 
ature [48,49]. 

The simulations in this work use a setup of two in- 
finite planes separated by the distance 2d. We call the 
direction between the two planes x and if not stated 
otherwise 2d is set to 64 lattice sites. In y direction 
we apply periodic boundary conditions. Here, 8 lattice 
sites are sufficient to avoid finite size effects since there 
is no propagation in this direction, z is the direction of 
the flow with our channels being 512 lattice sites long. 
At the beginning of the simulation (t — 0) the fluid is 
at rest. We then apply a pressure gradient VP in the z- 
direction to generate a planar Poiscuillc flow. Assuming 
Navier's boundary condition, the slip length (3 is mea- 
sured by fitting the theoretical velocity profile as given 
by equation [2]in flow direction (v z ) at position x, to the 
simulated data via the slip length /3. We validate this 
approach by comparing the measured mass flow rate 
/ T]v(x)dx to the theoretical mass flow without bound- 
ary slip and find a very good agreement. The dynamic 
viscosity \i as well as the pressure gradient ^ needed to 
fit equation ^ are obtained from our simulation data. 

In [47] , we show that this model creates a larger slip 
P with stronger interaction, namely larger g Q!Wa ii and 
larger i] wa ". The maximum available slip length mea- 
sured is 5.0 in lattice units. For stronger repulsive po- 
tentials, the density gradient at the fluid- wall interface 
becomes too large, causing the simulation to become 
unstable. At lower interactions the method is very sta- 
ble and the slip length f3 is independent of the distance 
d between the two plates and therefore independent of 
the resolution. We also show that the slip decreases 
with increasing pressure since the relative strength of 
the repulsive potential compared to the bulk pressure is 
weaker at high pressure. Therefore, the pressure reduc- 
tion near the wall is less in the high pressure case than 
in the low pressure one. Furthermore, we demonstrate 
that j3 can be fitted with a semi analytical model based 
on a two viscosity model. 
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Fig. 1 Slip length (3 versus bulk velocity v for different fluid- 
wall interactions r? wa11 . (3 is independent of v and only depends 
on r/ wa11 [47]. All units are expressed in lattice units throughout 
this paper if not stated otherwise. 



3 Roughness induced apparent slip 

If typical length scales of the experimental system are 
comparable to the scale of surface roughness, the ef- 
fect of roughness cannot be neglected anymore. Fig- 
ure[2|left) shows a typical example of a simulation setup: 
Poiseuille flow between two rough surfaces. The surface 
is generated using a random number generator to ran- 
domly choose the height of the obstacles at every dis- 
crete surface position. As can be observed in the figure, 
the stream lines of the flow are getting disturbed or 
trapped between the obstacles at the surfaces. In this 
section we show that an apparent boundary slip can 
have its origin in the misleading assumption of perfectly 
smooth boundaries. 



We study the dependence of the slip length /3 on 
the flow velocity for a wide range of velocities of more 
than three decades as it can be seen in Fig. [I] and 
in [47]. In the figure, we show data for different fluid- 
wall interactions < rj wal1 < 2.0 and flow velocities 
from 10~ 4 < v < 10 _1 . For simplicity we restrict our- 
selves to .g Qi waii = 0.08 which is a suitable value found 
from parameter studies given in [47]. Within this re- 
gion we confirm the findings of many steady state ex- 
periments [6], i.e., that the slip length is independent of 
the flow velocity and only depends on the wettability 
of the channel walls. Some dynamic experiments, how- 
ever, find a shear rate dependent slip [21,50]. These 
experiments often utilize a modified AFM as described 
in the introduction to detect boundary slippage. Since 
the slip length is found to be constant in our simula- 
tions after sufficiently long simulation times, we can- 
not confirm these results. However, it has been pro- 
posed by various authors that this velocity dependence 
is due to non-controlled effects such as impurities or 
surface nanobubbles. In simulations we can only find a 
shear rate dependence if the system has not yet reached 
the steady state or if time-dependent accelerations are 
present [51]. 

Our mesoscopic approach is able to reach the small 
flow velocities of known experiments and reproduces 
results from experiments and other computer simula- 
tions, namely an increase of the slip with increasing 
liquid-solid interactions, the slip being independent of 
the flow velocity, and a decreasing slip with increasing 
bulk pressure. In addition, within our model we develop 
a semi-analytic approximation of the dependence of the 
slip on the bulk pressure as described in [47]. 




Fig. 2 Left: a typical simulated system. Poiseuille flow between 
two rough surfaces showing random surface variations. Stream- 
lines depict a two dimensional cut and illustrate the parabolic 
velocity profile. This profile is distorted in the vicinity of the 
rough surfaces [52]. Right: the effective boundary height h e g is 
found between the deepest valley at ft. m j n and the highest peak 
at Vax- It corresponds to an effective channel width d e ff- Ra de- 
notes the average roughness and the maximum distance between 
the plates e£ max is kept constant [52], 



The influence of surface variations on the slip length 
[3 has been investigated by numerous authors. It was 
demonstrated by Richardson that roughness leads to 
higher drag forces and thus to no-slip on macroscopic 
scales. He has shown that if on a rough surface even 
a full-slip boundary condition is applied, one obtains 
a flow speed reduction near the boundary resulting in 
a macroscopic no-slip boundary condition [53]. An ex- 
perimental confirmation was later presented by McHale 
and Newton [54]. Molecular dynamics (MD) simula- 
tions of Couette flow between sinusoidal walls have been 
presented by Jabbarzadeh et al. [55]. They found that 
slip appears for roughness amplitudes smaller than the 
molecular length scale [55]. Sbragaglia et al. applied 
the LB method to simulate fluids in the vicinity of mi- 
crostructured hydrophobic surfaces [56], Al-Zoubi et al. 
demonstrated that the LB method is well applicable 
to reproduce known flow patterns in sinusoidal chan- 
nels [57], and Varnik et al. [58,59] have shown that 
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even in small geometries rough channel surfaces can 
cause flow to become turbulent. 

Recently, we presented the idea of an effective wall 
for rough channel surfaces [60]. Here, we investigate the 
influence of different types of roughness on the position 
of the effective boundary. Further, we show how the 
effective boundary depends on the distribution of the 
roughness elements and how roughness and hydropho- 
bicity interact with each other [52]. Lecoq and cowork- 
ers [61] performed experiments with well defined rough- 
ness, and developed a theory to predict the position of 
the effective boundary. In the experiments they utilised 
a laser interferometer to measure the trajectory of a 
colloidal sphere, and thereby determined the lubrica- 
tion force and an effective boundary position. The used 
geometry consists of grooves with a triangular profile. 
For a theoretical description the boundary is expressed 
in a Fourier series that gives the boundary condition for 
the Laplace equation. From this an effective boundary 
can be derived by a fast converting series. 

In this paper, we revise our previous achievements 
and compare them with the theoretical and experimen- 
tal results of Lecoq and coworkers [61]. 

Again, Poiseuille flow measurements are utilized to 
investigate the effect of interest. The rough surfaces are 
characterised by the highest point of one plane (/i ma x), 
the position of the deepest valley (/i m i n ) and the arith- 
metic average of all surface heights giving the average 
roughness R a . In the case of symmetrical distributions 

We get R a = ft-max/2. 

The position of the effective boundary h e g can be 
found by fitting the parabolic flow profile via the dis- 
tance d e s. With P set to we obtain the no-slip case. 
To obtain an average value for the effective distance 
between the planes d G ff , a sufficient number of individ- 
ual profiles at different positions z are taken into ac- 
count. The so found d e g gives the position of the effec- 
tive boundary and the effective height h c g of the rough 
surface is then defined by <f max — d c s (see Fig. [2] left). 

We show that the position of the effective boundary 
height is depending on the shape of the roughness ele- 
ments, i.e., for strong surface distortions it is between 
1.69 and 1.90 times the average height of the roughness 
Ra = ^max/2 [60]. In Fig [3] we plot the effective bound- 
ary positions of different geometries, i.e. randomly dis- 
tributed grooves with a square profile and grooves with 
a triangular profile. The results for the triangular ones 
match with the theoretical value of Lecoq et al. [61] for 
a similar geometry. 

By adding an additional distance between roughness 
elements, h c g decreases slowly, so that the maximum 
height is still the leading parameter. We are also able 
to simulate flow over surfaces generated from AFM data 



CD 

_cz 

CD 
> 

O 
CD 



18 
16 
14 
12 
10 h 

8 

6 

4 

2 





triangles * 

blocks ■ 

random ♦ 
Lecoq et al 




1 2 3 4 5 6 7 8 
average roughness R a 



Fig. 3 Simulated effective height h e g versus R a for different 
surface geometries. The triangular shape matches the theoretical 
results of Lecoq et al. [61] for a similar geometry. 



of gold coated glass used in microflow experiments by 
O.I. Vinogradova and G.E. Yakubov [62]. We find that 
the height distribution of such a surface is Gaussian 
and that a randomly arranged surface with a similar 
distribution gives the same result for the position of 
the effective boundary although in this case the heights 
are not correlated. 
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Fig. 4 Simulated effective height h e g versus R a for gold coated 
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We can tune the width of the distribution a and 
the average height R a . By scaling a with R a we ob- 
tain geometrically similar geometries. This similarity is 
important because the effective height h c g scales with 
the average roughness in the case of geometrical simi- 
larity [60]. We investigate Gaussian distributed heights 
with different widths a and find that the height of the 
effective wall depends linearly on a in the observed 
range [52]. Further, we find that the slip diverges as 
the amplitude of the roughness increases and the flow 
field gets more restricted which highlights the impor- 
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tance of a proper treatment of surface variations in very 
confined geometries [60]. 

4 Structured surfaces with entrapped 
microbubbles 

A natural continuation of our previous works on rough- 
ness induced apparent boundary slip and the collab- 
oration mentioned above is the analysis of flow along 
superhydrophobic surfaces [63] . While in typical exper- 
iments, slip lengths of a few tens of nanometers can be 
observed, it would be preferable for technical applica- 
tions to increase the throughput of fluid in a microchan- 
nel, i.e., to obtain substantially larger slip. Superhy- 
drophobic surfaces are promising in this context, since 
it has been recently predicted [64] and experimentally 
reported [65] that the so-called Fakir effect or Cassie 
state considerably amplifies boundary slippage. Using 
highly rough hydrophobic surfaces such a situation can 
be achieved. Instead of entering the area between the 
rough surface elements, the liquid remains at the top of 
the roughness and traps air in the interstices. Thus, a 
very small liquid-solid contact area is generated. 

Steinberger et al. utilized surfaces patterned with a 
square array of cylindrical holes to demonstrate that 
gas bubbles present in the holes may cause a reduced 
slip [66]. Numerically, they found even negative slip 
lengths for flow over such bubble mattresses, i.e., the 
effective no-slip plane is inside the channel and the bub- 
bles increase the flow resistance. In this section we con- 
sider negative slip lengths on bubble surfaces and also 
discuss the question of shear-rate dependent slip. In 
particular we show that microbubbles can generate a 
shear-rate dependence. 

Our simulations utilize the single component multi- 
phase LB model by Shan and Chen [46] which enables 
simulations of liquid-vapor systems with surface ten- 
sion. We are not aware of further lattice Boltzmann 
simulations to study the flow over a bubble mattress. 
However, a number of authors has applied various LB 
multiphase and multicomponent models to study the 
properties of droplets on chemically patterned and su- 
perhydrophobic surfaces [67-70] . The flow in our system 
is confined between two parallel walls. One of the walls 
is patterned with holes and vapor bubbles are trapped 
to these holes. The other wall is smooth and moved 
with velocity uo. Steinberger et al. [66] presented finite- 
element simulations of flow over rigid "bubbles" by ap- 
plying slip boundaries at static bubble surfaces. The 
LB method allows the bubbles to deform if the viscous 
forces are high enough compared to the surface tension. 
We are also interested in how surface patterning affects 
the slip properties of these surfaces, and how bubbles 



could be utilized to develop surfaces with special prop- 
erties for microfluidic applications [63]. 




Fig. 5 A visualisation of the simulation setup (left): the lower 
surface is patterned with holes, while the upper surface is moved 
with velocity uo- Right: the slip length /3 as a function of pro- 
trusion angle tp. A unit cell of each array is shown in insets and 
corresponding results are given by triangles (rhombic array), di- 
amonds (rectangular array), and circles (square array). The inset 
in the top- left corner shows the definition of ip [63] . 

The distance between walls is d = 1 /im (40 lattice 
nodes) in all simulations, and the area fraction of holes 
is 0.43. A unit cell of the regular array is included in 
a simulation and periodic boundary condition are ap- 
plied at domain boundaries. The bubbles are trapped 
to holes by using different wettabilities for boundaries 
in contact with the main channel and with the hole. 
The protrusion angle ip (see Fig.[5]for definition) is var- 
ied by changing the liquid's bulk pressure. The effective 
slip length is (3 = fiuo/a — d, where a — /idv/dz is the 
shear stress acting on the upper wall and /i the dynamic 
viscosity of the liquid. 

We investigate the effect of a modified protrusion 
angle and different surface patterns by using square, 
rectangular, and rhombic bubble arrays. The cylindri- 
cal holes have a radius a = 500 nm and the area fraction 
of the holes is equal in all cases. The shear rate is such 
that the Capillary number Ca = [iaG s /j = 0.16. Here, 
G s and 7 are the shear rate and surface tension, respec- 
tively. A snapshot of a simulation is shown in the left 
part of Fig. [5] and the slip lengths obtained are shown 
in the right part. The observed behavior is similar to 
that reported in [66] , where a square array of holes was 
studied. In particular, we observe that when (p is large 
enough (3 becomes negative. Moreover, when the pro- 
trusion angle equals zero, the slip length is maximised 
and the highest possible throughput in a microchannel 
is obtained. The behavior of the slip length can be ex- 
plained by thinking of an increased surface roughness 
if the protrusion angle is larger or smaller than zero. 
Since the area fraction of the bubbles is the same in all 
three cases, our results clearly indicate that slip proper- 
ties of the surface can be tailored not only by changing 
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the protrusion angle but also by the array geometry. 
In the presented study, the highest slip lengths are ob- 
tained for the rhombic unit cell and it is current work 
of progress to investigate the influence of the array ge- 
ometry in more detail. Recently, our findings have been 
confirmed theoretically by Davis and Lauga [71]. 

Next, the shear-rate dependence of the slip length 
is investigated. As the shear rate and thus the viscous 
stresses grow the bubbles are deformed (see Fig. [6] left) 
and the flow field is modified. In the central part of 
Fig. [6] we show the simulated slip length as a func- 
tion of the Capillary number for three different pro- 
trusion angles. The Capillary numbers chosen are in 
higher end of the experimentally available range. Our 
results show shear-rate dependent slip, but the behavior 
is opposite to that found in some experiments: in fact, 
the slip lengths measured by us decrease with increas- 
ing shear due to a deformation of the bubbles. In the 
experiments, surface force apparatuses are used (see, 
e.g., Ref. [21]), where a strong increase in the slip is ob- 
served after some critical shear rate. This shear-rate de- 
pendence has been explained, e.g., with formation and 
growth of bubbles [12,72]. In our simulations, there is 
no formation or growth of the bubbles as we only sim- 
ulate a steady case for given bubbles. The experiments 
on the other hand are dynamic. However, our results 
indicate that the changes in the flow field which occur 
due to the deformation of the bubbles cannot be an 
explanation for the shear-rate dependence observed in 
some experiments. Our results are consistent with [60] 
and the previous section, where it is shown that smaller 
roughness leads to smaller values of a detected slip. In 
the present case, the shear reduces the average height of 
the bubbles and thus the average scale of the roughness 
decreases as well. 

Finally, we consider a surface patterned with grooves. 
Cylindrical bubbles protrude to the flow channel from 
these holes with protrusion angle <p — 72°, and the 
area fraction of slots is 0.53. We apply shear both par- 
allel and perpendicular to the slots. The slip length is 
strongly dependent on the flow direction [63]. For par- 
allel flow the slip length is positive, but for the perpen- 
dicular case it becomes negative. Flow direction affects 
also greatly on the shear-rate dependence (cf. Fig. [6] 
right). When flow is parallel to the grooves no shear- 
rate dependence is observed, but for the perpendicular 
case this dependence is similar to that seen on hole ar- 
rays. These results can be understood on the basis of 
deforming bubbles. For perpendicular flow the bubbles 
are able to deform, but for the parallel case the bubbles 
retain their shape regardless of the shear rate. 



5 Conclusion 

In this paper we review applications of the lattice Boltz- 
mann method to microfluidic problems. The main focus 
of the paper is on our own research related to the vali- 
dation of the no-slip boundary condition. By introduc- 
ing a model for hydrophobic fluid-surface interactions 
and studying pressure-driven flow in microchannels, we 
show that an experimentally detected slip can have its 
origin in hydrophobic interactions, but is constant with 
varied shear rates and decreases with increasing pres- 
sure. Another effect that was not fully understood so 
far is the influence of surfaces roughness. We are able 
to apply our simulations to surface data obtained from 
AFM measurements of experimental samples. We show 
that ignoring roughness can lead to large errors in a 
detected slip. In fact, we propose that roughness alone 
could often be the reason for apparent boundary slip. 

Microscalc bubbles at surfaces allow to tailor the 
slip properties of a surface. Such a surface with bubbles 
may yield negative slip, i.e., increased resistance to flow, 
if bubbles are strongly protruding to the channel. The 
lattice Boltzmann simulations capture the deformabil- 
ity of bubbles and thus allow to study the influence 
of the shear rate on the deformation of the interface 
and it's effect on the measured slip. We find that the 
slip decreases with increasing shear rate demonstrating 
that shear induced bubble deformation cannot explain 
recent experimental findings where slip increases with 
increasing shear rate. 

In the current review, we also demonstrate the suit- 
ability of the lattice Boltzmann method for modeling 
microfluidic applications: in contrast to molecular dy- 
namics, it is able to reach experimentally available time 
and length scales. This allows one to compare simula- 
tion results to experimental data directly as demon- 
strated in the case of simulations of flow along surface 
data obtained from AFM measurements of "real" sam- 
ples. 
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Fig. 6 The left figure shows a snapshot of a bubble deformed by shear flow. In the centre, the slip length as a function of the capillary 
number for a square array of bubbles with three different protrusion angles, ip = 63°, 68°, and 71° (from uppermost to lowermost) is 
shown. The inset shows cross sections of liquid-gas interfaces for four capillary numbers [63] . The right figure shows the slip length as 
a function of capillary number for a surface with cylindrical bubbles. Circles denotes the values for flow parallel to the bubbles and 
diamonds for the perpendicular direction. 



